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Abstract 

Using non-equilibrium molecular dynamics method(NEMD), we have found that the thermal 
conductivity of multilayer graphene nanoribbons monotonously decreases with the increase of the 
number of layers, such behavior can be attributed to the phonon resonance effect of out-of-plane 
phonon modes. The reduction of thermal conductivity is found to be proportional to the layer size, 
which is caused by the increase of phonon resonance. Our results are in agreement with recent 
experiment on dimensional evolution of thermal conductivity in few layer graphene. 
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The fresh comings of carbon family, the two-dimensional (2D) graphene and quasi-one- 
dimensional (ID) graphene nanoribbon (GNR), have attracted strong interest due to their 
fundamental physical properties and potential applications in nanoelectronic devices{l|- 
1 Experimental developments have enabled the growth of the h,gh quality multilayer 
graphenep] and the control of GNR edge geometries|6|. As we know, the thermal property 
is a crucial aspect that determines the application of a material in the nanoelectronics[7|. In 
recent years, the thermal properties of monolayer graphene and GNRs have aroused much 
attention 8-1^. Very recently, the thermal conductivity of few layer graphene has been also 



studied experimentally and theoretically 



5fl. However, the thermal conductivity of multilayer 



GNRs have not been well investigated so far. Its clarification is beco ming desirable due to 



the forthcoming application of multilayer GNRs in nanoelectronics 



18|-|20|. 



On the other hand, the mechanism of thermal conduction on the nanoscale is currently 
a controversial issue^, ^2^. Investigating the dimensional evolution of thermal conductivity 
could provide a new insight to clarify the fundamental mechanism on nanoscale. The thermal 
conductivity evolution from 2D to 3D had been investigated experimentally jl]], while the 
evolution from ID to higher dimensions is still not well studied. The multilayer GNRs 
are ideal systems for such investigation, whose dimensional evolution can be easily realized 
through the layer number and size variation. 

In this letter, we employ the non-equilibrium molecular dynamics (NEMD) method to 
investigate layer and size effect on thermal conductivity of multilayer GNRs with different 
edge shapes, such as the armchair multilayer GNRs and the zigzag multilayer GNRs. The 
number of layers varies from 1 to 4 and the width of layers varies from 1 to 10 nm. We 
find that both the number of layers and size have strong influence on thermal conductivity 
of multilayer GNRs. The intrinsic mechanism of thermal conductivity variation is further 
explored through a phonon resonance model. 

The structure of multilayer GNRs are constructed based on the theoretical prediction of 



bilayer GNRs which has a small deviation from Bernal stacking 



tion, Tersoff potential 



22|. In the NEMD simula- 



23] is utilized to describe the in-plane C-C bonding interactions and 



Lennard- Jones potential is used to describe the intra-plane van der Waals interactions 



2J. 



We use the velocity Verlet method to integrate equations of motion with a fixed time step of 



1 fs 
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25M27]. On each layer of multilayer GNR, fixed boundary condition is_ implemented 

17|. Next to 



with the atoms at the left and right ends fixed at their equilibrium positions 



16 
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the boundaries, the adjacent two cells of atoms are coupled to Nose — Hoover thermostats 
with temperatures 310 K and 290 K, respectively. The thermal conductivity for layer i can 
be calculated directly from the well known Fourier law Ki = Jidf(ATwh), where AT =20 
K is the temperature difference between two thermostats, Jj is the heat flux from the heat 
bath to the system , which can be obtained via calculating the power of heat baths|25|, d is 
the length, w is the width, h ( 0.144 nm ) \v\ is the thickness of a monolayer GNR. All the 
calculated thermal conductivities are obtained by averaging about 3 ns after 2 ns to estab- 
lish a stable temperature gradient along the length direction. Thus the thermal conductivity 
of multilayer GNR can be defined as k = ^2 Ki/n with n being the number of layers. In 
addition, all the GNR structures have been optimized before NEMD simulation. We also 
define a multilayer GNR with N carbon-chains in width to be represented as iV-armchair 
multilayer GNR or iV-zigzag multilayer GNR, depending on the specific edge shapes jl7J. 

We first investigate the number of layers dependence of thermal conductivity of 20- 
armchair multilayer GNR and 10-zigzag multilayer GNR. Similar to that of monolayer 
GNRs, the 20- armchair multilayer GNRs have much lower thermal conductivity than 10- 
zigzag multilayer GNRs (see Fig. 1), implying a universal edge-shape dependence of thermal 
conductivity in the GNR family. Moreover, with the number of layers increasing, thermal 
conductivity of both the armchair multilayer GNR and zigzag multilayer GNR monotonously 
decreases. When the number of layers gets to 4, the thermal conductivity is reduced to 123 
Wm" 1 K _1 and 308 Wm _1 K _1 for 20-armchair multilayer GNR and 10-zigzag multilayer 
GNR, respectively. Comparing with that of monolayer 20-armchair GNR (195 Wm^K" 1 ) 
and 10-zigzag GNR (495 Wm~ 1 K~ 1 ), the reduction of thermal conductivity of armchair 
multilayer GNR and zigzag multilayer GNR gets to 40%, showing an obvious dependence 
on number of layers. This indicates that the crossplane coupling is enhanced with number 
of layers increasing, and it plays an important role in the evolution of thermal conductivity 
from monolayer GNR to multilayer GNR. This is in agreement with recent experiment in few 
layer graphene, where a 67% reduction of thermal conductivity is observed as the number 
of atomic plane increasing from 1 to 



These results can be further understood by the coupling mechanism |28|], there exists a 
competitive mechanism on thermal conductivity in a coupling system: the phonon-resonance 
effect that decreasing thermal conductivity and phonon-band-up-shift effect that increasing 



thermal conductivity 



28[. The strength of phonon resonance can be described by the res- 
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onance angle which is determined by atomic mass and spring constant of two coupled 
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systems. When \I/ is small(\l/ < — ), the variation of thermal conductivity is dominated by 
the phonon-band-up-shift effect; when \1> comes to large(\I / > — ), the thermal conductivity 
is dominated by the phonon resonance effect. The thermal conductivity reduction of mul- 
tilayer GNRs with the increase of number of layers can be accounted for this mechanism. 
For the multilayer GNRs, the atomic mass and spring constant of the coupling layers are 

7T 

totally equivalent, thus \l/ = — and the phonon-resonance effect plays a dominant role on 
thermal conductivity variation. Therefore, the thermal conductivity of multilayer GNRs 
monotonously decreases with the number of layers increasing which induces more and more 
intensive phonon resonance. 

In order to identify the phonon-resonance effect in multilayer GNRs, we freeze the out- 
of-plane atomic vibration of multilayer GNRs and re-calculate their thermal conductivity. 
For simplicity, we consider the bilayer GNRs. Here we only present the calculated thermal 
conductivity values for the zigzag bilayer GNR, the results for armchair bilayer GNR are 
qualitatively similar. Two cases are considered: out-of-plane vibration of the top layer is 
frozen(constraint 1) and out-of-plane vibration of both layers is frozen(constraint 2). 

As shown in Table 1, freezing out-of-plane atomic vibration would considerably change 
the layer's thermal conductivity. If only top layer's out-of-plane vibration is frozen, thermal 
conductivity of top layer would increase by 40% (from 334 Wm~ 1 K~ 1 to 467 Wm^K -1 ) 
while thermal conductivity of bottom layer is nearly unaffected by the artificial constraint. 
If the out-of-plane vibration is frozen in both layers, thermal conductivity of bilayer zigzag 
GNR almost equals to that of monolayer zigzag GNR. This indicates that, in the multilayer 
GNRs, the phonon resonance of stacking layers is mainly from the coupling between the 



29| and the out-of-plane phonon modes that 



crossplane out-of-plane ZO' phonon modes 
propagate in the basal plane. 

For the purpose of investigating finite size effect, we calculate thermal conductivity of 
both monolayer and bilayer GNRs with various width of GNRs. Figure 2 shows the obtained 
thermal conductivity of monolayer and bilayer GNRs, whose width ranges from 1 nm to 10 
nm. As one can see, the thermal conductivity of zigzag bilayer GNR increases firstly and 
then turns to decrease with the width increasing, while the armchair bilayer GNR's thermal 
conductivity monotonously increases with the width increasing. The width dependence of 
thermal conductivity of bilayer GNRs is very similar to that of monolayer GNRs (Fig. 2). 
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This phenomena can be attributed to the competition between edge-localized phonon effect 



id- 



and the umklapp phonon scattering effect with the width of ribbon variation 

In addition, the difference of thermal conductivity between monolayer and bilayer GNRs 
also exhibits an obvious layer-size dependence. For the armchair GNRs (zigzag GNRs), the 
reduction of thermal conductivity monotonously increases from 31% (30%) to 35% (32%) 
with the width increasing from 1 nm to 10 nm. This means the difference of thermal 
conductivity between monolayer GNR and multilayer GNR is proportional to the layer 
size, because the phonon resonance strength between different layers is proportional to the 
number of total phonon modes which is in turn corresponding to the size of ribbon. Our 
results on thermal conducitivity reduction of multilayer GNRs is smaller than that of few 
layer graphenesQ. Above results indicate that the difference in dimension evolution of 
thermal conductivity between multilayer graphene and multilayer GNRs comes from the 
finite size effect. 

In conclusion, we have investigated the thermal conductivity of multilayer GNRs using 
the NEMD method. Comparing to monolayer GNR, the thermal conductivity of multilayer 
GNRs has a significant reduction which is in agreement with available experiments p]. Based 
on the phonon resonance model, we propose that the reduction of thermal conductivity is 
attributed to the resonance of out-of-plane phonon modes. Moreover, the difference between 
thermal conductivities of GNRs with different number of layers is found to be proportional 
to the layer size, which is directly determined by the number of the total phonon modes. 
The present studies suggest that the thermal conductivity of multilayer GNRs can be ma- 
nipulated by changing the number and size of layers, which provides potential applications 
of multilayer GNR-based materials in future nanodevices. 
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puter Center of Fudan University. 
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Table I: Thermal conductivity k of bilayer zigzag GNRs with different constraints. The length and 
width of the monolayer and bilayer zigzag GNR is about 10 nm and 3 nm. 



Type of GNR 

Bilayer zigzag GNR 
Bilayer zigzag GNR 
Bilayer zigzag GNR 
Monolayer zigzag GNR 



Constraint 

Free 
Constraint 1 
Constraint 2 

Free 



Thermal Conductivity (Wm K ] 
Top Layer Bottom layer 

323 335 
330 467 
469 480 
490 



Constraint 1: Out-of-plane modes of top layer are frozen Constraint 2: Out-of-plane modes of both layers are frozen 
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Figure 1: Thermal conductivity n as a function of the number of layers with length and width of 
5 nm and 2 nm. k decreases monotonously with number of GNR layers increasing, indicating the 
enhancement of intra-layer phonon coupling. 
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Figure 2: Thermal conductivity n as a function of the width of monolayer or bilayer armchair GNR 
and zigzag GNR with fixed length of 10 nm. The behavior of k for monolayer and bilayer is similar 
and the difference of K between monolayer and bilayer GNRs increases with the width of GNRs 
increasing. 
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